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Abstract 

In this work a two-particle irreducible (2PI) closed-time-path (CTP) effective action is used to describe 
the nonequilibrium dynamics of a Bose Einstein condensate (BEC) selectively loaded into every third site 
of a one-dimensional optical lattice. The motivation of this work is the recent experimental realization of 
this system at National Institute of Standards and Technology (NIST) where the placement of atoms in an 
optical lattice is controlled by using an intermediate superlattice. This patterned loading method is a useful 
technique for the proposed implementation of lattice-based AMO quantum computing. This system also 
serves to illustrate many basic issues in nonequilibrium quantum field theory pertaining to the dynamics 
of quantum correlations and fluctuations which goes beyond the capability of a mean field theory. By 
numerically evolving in time the initial state configuration using the Bose- Hubbard Hamiltonian an exact 
quantum solution is available for this system in the case of few atoms and wells. One can also use it to test out 
the various approximation methods constructed. Under the 2PI CTP scheme with this initial configuration, 
three different approximations are considered: a) the Hartree-Fock-Bogoliubov (HFB) approximation, b) 
the next-to- leading order l/M expansion of the 2PI effective action up to second order in the interaction 
strength and c) a second order perturbative expansion in the interaction strength. We present detailed 
comparisons between these approximations and determine their range of validity by contrasting them with 
the exact many body solution for a moderate number of atoms and wells. As a general feature we observe 
that because the second order 2PI approximations include multi-particle scattering in a systematic way, 
they are able to capture damping effects exhibited in the exact solution that a mean field coUisionless 
approach fails to produce. While the second order approximations show a clear improvement over the HFB 
approximation our rmmerical result shows that they do not work so well at late times, when interaction 
effects are significant. 



1 Description of the Problem 

Bose-Einstein condensate (BEC) loaded into an optical lattice has provided an interesting arena for the study 
of quantum coherence and fluctuation phenomena in many body physics. Spectacular progress in experimental 
studies have been able to achieve regimes where standard mean field techniques used to describe weakly inter- 
acting atoms are not generally applicable. The description of the evolution of condensates far from equilibrium 
has also gained considerable importance in matter-wave physics, motivated by recent experimental realization 
of colliding and collapsing condensates. In this paper we want to investigate the dynamics of a Bose-Einstein 
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condensate at zero-temperature (T=0) loaded every three sites in a one dimensional optical lattice, recently 
realized by the NIST group 1 . 

In that experiment a combination of two independently controlled lattices was used to prepare the condensate 
into every third site of a single lattice. Briefly the procedure consists of loading a condensate into the ground 
band of a lattice with periodicity 3o, so that the condensate is well localized in the potential minima of this 
lattice. A second lattice of periodicity a which is parallel to the first lattice is then ramped up so that the 
superimposed light potentials form a super-lattice of period 3a. If both lattices are in phase, then the addition 
of the new lattice will not shift the locations of the potential minima from those of the first lattice alone and the 
condensate will remain localized at these positions. Finally, by removing the first period 3a lattice on a time 
scale long compared to band excitations, but short compared to the characteristic time of dynamics within the 
lattice, the condensate will be left in every third site of the period a lattice. 

A condensate loaded in this way is not an eigenstate of the final period a lattice, and the condensate will 
continue to evolve in the final system. Outside the strongly correlated regime, a mean field approach was 
expected to give a good description of the condensate dynamics '2'. However, by comparison with numerical 
simulations of the exact solution of the many-body Hamiltonian it was found that even in the case when 
the kinetic energy is comparable with the on-site interactions, inter-atomic collisions play a crucial role in 
determining the quantum dynamics of the system and therefore a mean field coUisionless approach is only 
accurate for short times. 

Thus, in order to model the correct quantum dynamics of the system it is necessary to include properly 
scattering processes among particles. This task, however, is not easy for this particular system because contrary 
to the three dimensional dilute gas case, where many body effects introduce only a small change to the two- 
particle scattering properties in vacuum, the presence of the lattice and the low dimensionality of the system 
make the problem much less straightforward. 

To date most theoretical descriptions of nonequilibrium dynamics are based on the time dependent Gross- 
Pitaevskii equation coupled with extended kinetic theories that describe excitations near thermal equilibrium 
situations. These approaches usually rely on the contact-potential approximation of the binary interactions 
which restricts two-body collisions to low energy. However, system dynamics of most interest falls outside of 
this restricted domain; its description thus requiring new methods. Here to describe the far-from-equilibrium 
dynamics we adopt a OTP (Closed Time Path) ^ functional-integral formalism together with a 2PI (two- 
particle irreducible) ^ effective action approach to derive the equations of motion. This method has been 
generalized for and applied to the establishment of a quantum kinetic field theory [S] [HJ [7] with applications to 
problems in gravitation and cosmology [U IHj , particles and fields ^ , BEC El and condensed matter 
systems ^ '^^^^ ^ addressing the issues of thermalization and quantum phase transitions ^1 El El ■ 

In Section 2 we summarize the mean field results obtained in previous studies 0. In Section 3 we introduce 
the 2PI generating functional to construct the 2PI effective action r2 and Green's functions. In Section 4 we 
perform perturbative expansions on r2 and define the various approximation schemes. In Section 5 we introduce 
the CTP formalism. We then derive the equation of motion and discuss the results under each approximation 
scheme, starting with HFB in Section 6, followed by second order expansions in Section 7, which includes the 
1/Af expansion and the full second order expansion. The numerical implementation is discussed in Section 8. 
In Section 9 we present our results and determine the range of validity of the approximations by comparisons 
with the exact (numerical) solution. We end with the conclusions. 

2 Mean Field Theory 

The dynamics of an ultracold bosonic gas in an optical lattice can be described by a Bose-Hubbard model where 
the system parameters are controlled by laser light. For a one dimensional lattice the starting Hamiltonian is: 
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where and $J are the annihilation and creation operators at the site i which obey the canonical commutation 
relations for bosons. Here, the parameter U denotes the strength of the on-site repulsion of two atoms on the 
site i; the parameter denotes the energy offset of each lattice site due to an additional slow varying external 
potential such as a magnetic trap and J denotes the hopping rate between adjacent sites. Because the next-to- 
nearest neighbor amplitudes are typically two orders of magnitude smaller, tunneling to them can be neglected. 
The Bose -Hubard Hamiltonian should be an appropiate model when the loading process produces atoms in 
the lowest vibrational state of each well, with a chemical potential smaller than the distance of the first 
vibrationally excited state. This is the case of the experiment we want to study T. 

Here we consider only a one-dimensional homogeneous lattice with periodic boundary conditions. The 
transverse degrees of freedom and the harmonic confinement in the direction of the lattice are relegated to a 
future work. 

In the simplest version of the mean-field theory, we assume that all atoms share the same condensate wave 
function and therefore it is possible to replace the operator on the lattice site i by a c-number 4>i{t). Under 
this approximation corresponds to the time dependent mean value amplitude of the i*^ well. One might 
expect that this zero temperature mean field approach could give a good description of the dynamics of the 
atoms in the limit that quantum fluctuations can be completely neglected, that is, the regime when the ratio U 
to J is sufficiently small and the tunneling overwhelms the repulsion. 

By making the mean field ansatz in the Bose- Hubbard Hamiltonian it is possible to show that the 
amplitudes satisfy the Discrete Non Linear Schrodinger Equation (DNLSE), which in the case of zero external 
potential has the form: 

where we have assumed that the ^^(t) are normalized to one, that N is the total number of atoms and defined 

a new time scaled t = H-. 

h 

We treat a model case in which the initial occupancies of each third site are the same, and in which the 
condensate initially has a uniform phase. Thus at t = 0, the amplitudes ^^(t) are given by '^'3^(0) = 
'^'si+ii^) = '^'31+2(0) = 0; where / is the total number of lattice sites. For an infinite lattice, or one with 
periodic boundary conditions, the amplitudes for all initially occupied sites (j)3i evolve identically in time, and 
the amplitudes for the initially unoccupied sites satisfy 4>3i+ii''') = (t>3i+2i''') fo^' This allows us to reduce 

the full set of equations (O to a set of two coupled equations for (I)q{t) and 0i(t). 

The solutions |</>o(''")| and |'/>i(t)| are oscillatory functions whose amplitudes and common period, 2^(7), are 
determined by the parameter 7 = NUI/{3J). It is useful to qualitatively divide the dynamical behavior into 
two regimes 

The tunneling dominated regime (7 < 1): In this regime we find that the oscillation period is essentially 
constant, the role of interactions is relatively small, and the equations of motion are equivalent to those of a 
two-state Rabi problem. This system will undergo Rabi oscillations whereby atoms periodically tunnel from 
the initially occupied site into the two neighboring sites. For 7 = the period of oscillation is ^ . 

Interaction dominated regime: The effect of interactions on the mean-field dynamics is to cause the 
energies of the initially occupied sites to shift relative to those of the unoccupied sites. As 7 increases the 
tunneling between sites occurs at a higher frequency, but with reduced amplitude. The population of the 
initially occupied sites becomes effectively self-trapped by the purely repulsive pair interaction. 

To check the validity of the mean field approximation, we made comparisons with the exact many body 
solution for 6 atoms and three wells. We use a modest number of atoms and lattice sites for the comparisons 
due to the fact that the Hilbert space needed for the calculations increases rapidly with the number of atoms 

and wells. The exact solution was obtained by evolving an initial state (e"^/^ e^*° |0))(8) |0) (g) |0) with the 
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Bose-Hubbard Hamiltonian. The initial state represents just a coherent state with an average of N atoms in 
the initial populated well. (See Sec. 8.1 ). 

In Fig. 1 we plot the average population per well (^^l {t)^i {t)j and the condensate population per well 

(t)^ and compare them with the mean field predictions, i.e. for three different values of 7. The 

salient features observed in these comparisons are: 



Weakly interacting regime (7 — 0.2): The DNLSE gives a good description of the early time dynamics. 
The total population per well predicted by the mean field solution agrees with the exact solution. This is 
because the total condensate population remains big for the time in consideration. 



Intermediate regime (7 = 2): Quantum fluctuations lead to a non-trivial modulation of the classical 
oscillations. In this regime the ratio between interaction and kinetic energy is small enough to allow the atoms 
to tunnel but not too small to make interaction effects negligible, the mean field results are accurate only for a 
short time. The exact solution exhibits for this regime damped oscillations of the atomic population per well 
that drives the system towards equilibrium. There is also a very rapid reduction of the condensate population. 
These characteristic damping effects are caused by the different phase evolution of the quantum states that 
span the configuration space and collisions among atoms. Quantum scattering effects are crucial, even for rather 
early times. 



Strongly correlated regime (7 = 12): the system exhibits macroscopic quantum self-trapping of the 
population. Qualitatively, both the mean field and the exact solutions, agree in the sense that atoms are 
in both self-trapped in the initially populated wells due to interactions. However, the fast decrease of the 
condensate population and its subsequent revivals shown in the exact solutions give us an idea of the importance 
of correlation effects beyond mean field in this interaction dominated regime. Due to the fact that the evolution 
of the initially unoccupied wells is almost frozen by interactions, the collapse and revival of the condensate 
population can be mainly described by the time evolution of the n-particlc number states that expand the 
coherent state of the initially populated well. ^ 

Because our main interest is the tunneling dynamics we will focus on the intermediate regime, where the 
ratio U/J is small but interaction effects are not negligible. In this regime a perturbative expansion around U/J 
still makes sense. 



3 2PI Effective Action r{(j), G) 

The first requirement for the study of nonequilibrium processes is a general initial- value formulation depicting 
the dynamics of interacting quantum fields. The CTP or Schwinger-Keldysh effective action formalism |3] serves 
this purpose. The second requirement is to describe the evolution of the correlation functions and the mean 
field on an equal footing. The two particle irreducible (2PI) formalism j4j where the correlation functions 
appear also as independent variables serves this purpose. By requiring the generalized (master) CTP effective 
action ^ to be stationary with respect to variations of the correlation functions an infinite set of coupled 
(Schwinger-Dyson) equations for the correlation functions is obtained which is a quantum analog of the BBGKY 
hierarchy. The 2PI effective action produces two such functions in this hierarchy. In this section we shall focus 
on the 2PI formalism, and then upgrade it to the CTP version in the next section. 

^The partial revival of the condensate population can be understood as the interference of states approximated by a quadratic 
spectrum: E„ U/2n{n — 1). At first, the different phase evolution of the n-particle states leads to a collapse of the condensate 
population. However at integer values of trev = {U/h)~^ the phase factors add to an integer value of 27r, leading to a revival of 
the initial state 1181 . The revival is not complete because J is not exactly zero. The collapse time tc depends on the variance of 
the initial atomic distribution tc ~ (treu/o") (See Ref. ll9l - i:^Ul V For the initial Poissonian distribution tc ~ {U\/~N /h)~^ . 
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An adequate description of the dynamics of atoms loaded into an optical lattice is provided by the Bose- 
Hubbard Hamiltonian (Q, whose classical action is given in terms of the complex fields $j and <&* by 



i 

where, as before, i denotes the lattice position, J is the hopping rate and U is the interaction strength. We 
limit the analysis to the case when no external potential is present and include only nearest neighbor hopping. 
To compactify our notation we introduce <i>°(a = 1, 2) defined by 

$, = $1, = (4) 

In terms of these fields the classical action takes the form 

sm = J dtY,^ha,<s>nt)hdM{t) + j^ab^^+i w^'w - ^{lab^trntityf, (5) 

i 

where Af is the number of fields, which is two in this case, and summation over repeated field indices a,b — (1, 2) 
is implied, hat and aab are matrices defined as 

^Qb = * ^ 0^ ) " ( 1 ) 

After second quantization the fields are promoted to operators. We denote the expectation value of the 
field operator or mean field by (j)i{t) and the expectation value of the composite field by Gfj{t, t'). Physically, 
is the condensate population and the composite fields determine the fluctuactions around the mean 

field. 



m) = mt)), (7) 
nGt^{t,t') = {Tcm)K{t')) - {m)) ■ (8) 

The brackets denote taking the expectation value respect to the density matrix and Tc denotes time ordering 
along a contour C in the complex plane. 

All correlation functions of the quantum theory can be obtained from the effective action r[(/), G], the two 
particle irreducible generating functional for Green's functions parametrized by (j\{t) and the composite field 
Gfj{t,t'). To get an expression for the effective action we first define the functional Z[J,K] W as 



where we have introduced the following index lowering convention 

Xa = CTabX\ (10) 

The functional integral © is a sum over classical histories of the field <i>° in the presence of the local source Jia 
and the non local source Jiijab- The coherent state measure is included in 13$. The addition of the two-particle 
source term is what characterizes the 2PI formalism. 
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We define r[0, G] as the double Legendre transform of T4/'[J. K] such that 



6W[3,K] 
SW[3,K] 

Expressing J and K in terms of and G yields 
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(11) 

(12) 



T[cj),G] = W[3,K] - I dtY,J^a{m{t) -\j dtdt'Y,^Ut)'f>''j{t')^^Mt,t') 

i ij 

-| I dtdt'Y,Gi!;{t,t')\^,,abM. 

ij 

From this equation the following identity can be derived: 



ST[(I>,G] 



= - i^.a{t)+ j dt'Y,{^^3ad{t.t'))<f>j{t')^ , 



(13) 



(14) 
(15) 



In order to get an expression for r[0, G] notice that by using Eq.lO for VF[J,K] and placing it in Ea. ll4|l 
for T[(j), G], it can be written as 



exp( -r[0,G] 



n J Z?<i>'^cxp I J (^Sm + J dU 3^a{t) mt) - CW] 



(16) 



-\ j dUdt', ($»(t)K,,a,(t,t')$?(i')-'/'- WK,,a,(t,t')0-(i')) - \TrGK)j 



= n/ ^^'^expjj (^Sm^J dU 



<5r[0,G] 



TrG 



6G 



where Tr means taking the trace. For simplicity we have denoted J dt^^ by J dti. Defining the fluctuation 
field, Lpf ^ ~ we have 



G; ^] = 5[0 + ^] - / dU ^-^^^nt) ~\j dUdt', 



(17) 
(18) 



By introducing the classical inverse propagator iD ^{(/)) given by 



6 



iD,,abit,t') = ^-^^ (19) 



= [Sijhabdt + J{5i+ij + Si^ij)aab)S{t - t') 
U 



the solution ol the lunctional integro-dilFerential equation p7|l can be expressed as 

r[0,G] = S[<p] + ^TrlnG'^ + ^TrD-^{(l))G + r2[(l),G] + const. (20) 

The quantity T2 [(j), G] is conveniently described in terms of the diagrams generated by the interaction terms 
in S[(j) + f] which are of cubic and higher orders in ip. 

S.nd^ + ip] - / dt, {^,,it)ip^{t)f '^Jdt, ^tma{tW■it)^^b{t). (21) 



4M J ' y^^^y-'^^y-" _^ 

It consists of all two-particle irreducible vacuum graphs (the diagrams representing these interactions do not 
become disconnected by cutting two propagator lines) in the theory with propagators set equal to G and vertices 
determined by the interaction terms in S[(l) + if] . 

Since physical processes correspond to vanishing sources J and K, the dynamical equations of motion for 
the mean field and the propagators are found by using the expression (|2U|I in equations H14|l and (|15|l , and 
setting the right hand side equal to zero. This procedure leads to the following equations: 

habnMlit) = -J(0.+iJi) + '/'.-iaW) + ;^(0.dW0fW + G,;?c(i, W+ (22) 

^AiGuadit, t) + G,,da{t, t))4{t) 



and 



G^lbit^') = D,jab{t,t')-' -^^MtA')^ (23) 

Sr2[(b,G] 
^SGf^Ht,t'y 



^^Mt^') ^ ^^TT^HTTk- (24) 



Equation H23|l can be rewritten as a partial differential equation suitable for initial value problems by 
convolution with G. This differential equation reads explicitly 

KndtG^''{t,t') = -J(G^^i^.(t,t') + Gfy(t,i')) + ^(0.<i(i)0?(O)G^''(i,i')+ (25) 
— <j>Ut)G'i^{t,t')^,,{t)+^ J dt';j:Ut,nGt,{t'',t')+z5^'6,,S{t-t'). 

The evolution of 0° and G"** is determined by Eqs. (|22|l and (|25|l once r2[(/), G] is specified. 
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4 Perturbative Expansion of r2 (0,(7) and Approximation Schemes 



The diagrammatic expansion of T2 is illustrated in Fig. 2, where two and three-loop vacuum diagrams are 
shown. The dots where four lines meet represent interaction vertices. The expression corresponding to each 
vacuum diagram should be multiplied by a factor (— where I is the number of solid lines and s the 
number of loops the diagram contains. 

The action F including the full diagrammatic series for F2 gives the full dynamics. It is of course not feasible 
to obtain an exact expression for F2 in a closed form. Various approximations for the full 2PI effective action 
can be obtained by truncating the diagrammatic expansion for F2. Which approximation is most appropriate 
depends on the physical problem under consideration: 

4.1 The standard approaches 

1. Bogoliubov (One- loop) Approximation: 

The simplest approximation consists of discarding F2 altogether. This yields the so called Bogoliubov or 
one-loop approximation whose limitations have been extensively documented in the literature ( j21l I22| 'l. 

2. Time-dependent Hartree-Fock-Bogoliubov (HFB) Approximation: 

A truncation of F2 retaining only the first order diagram in U, i.e., keeping only the double- bubble 
diagram, Fig. 2a, yields equations of motion of 4> and G which correspond to the time dependent Hartree- 
Fock-Bogoliubov (HFB) approximation. This approximation violates Goldstone's theorem, but conserves 
energy and particle number [I'M 124' '25"). The HFB equations can also be obtained by using cumulant 
expansions up to the second order f28 in which all cumulants containing three or four field operators are 
neglected. The HFB approximation neglects multiple scattering. It can be interpreted as an expansion in 
terms of Ut/J, (where t is the time of evolution) and is good for the description of short time dynamics 
or weak interaction strengths. It will be described in Sec. 6. 

4.2 Higher order expansions 

We make a few remarks on the general properties of higher order expansions and then specialize to two approx- 
imations. 

2PI and Ladder Diagrams Since the work of Beliaev's |^ and Popov [30] it is well known in the literature 
(see for example Ref. [121 that including higher order terms in a diagrammatic expansion corresponds to 
renormalizing the bare interaction potential to the four-point vertex, thus accounting for the repeated scattering 
of the bosons. In Ref. |S2 the authors have shown explicitly for a homogeneous Bose gas that taking into account 
the two-loop contribution to the 2PI effective action leads to diagrams topologically identical to those found by 
Beliaev but with the exact propagator instead of the one-loop propagator. In the dilute gas limit, where the 
inter-particle distance is large compared with the s-wave scattering length, the ladder diagrams give the largest 
contribution to the four-point vertex. Every rung in a ladder contributes to a factor proportional to Um (In 
the presence of the lattice m should be replaced by the effective mass m* ( m* ^ It? /{Jo?)) with a the lattice 
spacing). The ladder resummation results in an effective potential which is called the T- matrix. To lowest 
order in the diluteness parameter , the T-matrix in 3-D systems can be approximated by a constant proportional 
to the scattering length (pseudopotential approximation). However this approximation is only valid in the weak 
interaction limit and neglects all momentum dependence which appear in the problem as higher order terms. In 
that sense the 2PI effective action approach allows us to go beyond the weakly interacting limit in a systematic 
way and to treat collisions more accurately. 
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Nonlocal Dissipation and Non-Markovian Dynamics Higher order terms lead to nonlocal equations 
and dissipation. The presence of nonlocal terms in the equations of motion is a consequence of the fact that 
the 2PI effective action really corresponds to a further approximation of the master effective equation 6 . The 
2PI effective action is obtained by the slaving of the three point function C3 to the mean field and G with a 
particular choice of boundary conditions. See Ref. (HI for further details. 

Non-Markovian dynamics is a generic feature of the nPI formalism which yields integro-differential equations 
of motion. This makes numerical solution difficult, but is a necessary price to pay for a fuller account of the 
quantum dynamics. Many well acknowledged approaches to the quantum kinetics of such systems adopt either 
explicitly or implicitly (or at the end what amounts to) a Markovian approximation 34 . It assumes that 
only the current configuration of the system, but not its history, determines its future evolution. Markovian 
approximations are made if one assumes instantaneous interactions, or in the kinetic theory context that the 
time scales between the duration of binary collisions Tq and the inverse collision rate are well-separated. 
In the low kinetic energy, weak interacting regime the time between collisions (or mean free path) is long 
compared to the reaction time (or scattering length): Tc >> To- The long separation between collisions and 
the presence of intermediate weak fluctuations, allow for a rapid decay of the temporal and spatial correlations 
created between collision partners, which one can use to justify the Markovian approximation. However, in 
the problem at hand, owing to the presence of the lattice which confines the atoms to the bottom of the wells 
with enhanced interaction effects, the low dimensionality of the system and the far-from-equilibrium initial 
conditions, non-Markovian dynamics needs to be confronted squarely. That is the rationale for our adoption of 
the CTP 2PI scheme. Now, for the specifics: 

1. Second Order expansion: 

A truncation retaining diagrams of second order in U containing besides the double-bubble, also the setting- 
sun and the basket-ball Fig. 2. By including the setting-sun and the Basket-ball in the approximations 
we are taking into account two particle scattering processes IZ| ■ Second order terms lead to integro- 
differential equations which depend on the time history of the system. 

2. Large- A/" approximation 

The 1 /J\f expansion is a controlled non-perturbative approximation scheme which can be used to study 
non-equilibrium quantum field dynamics in the regime of strong interactions |lHI Il5| . In the large A/" 
approach the field is modeled by Af fields and the quantum field generating functional is expanded in 
powers of 1/Af. In this sense the method is a controlled expansion in a small parameter but unlike 
perturbation theory in the coupling constant, which corresponds to an expansion around the vacuum, the 
large Af expansion corresponds to an expansion of the theory about a strong quasiclassical field. 

In this work we will derive the equations of motion and perform numerical calculations up to the second 
order in the coupling constant U. This will enable us to determine the range of validity of three types of 
approximations described above, namely, a) the HFB, b) the full second order and c) NLO large Af expansion 
up to second order in U (in the figures we use the shorthand HFB, 2nd and 1/AA respectively) by comparison 
with the exact many body solution for a moderate number of atoms and wells. 

5 CTP formalism 

In order to describe the non-equilibrium dynamics we will now specify the contour of integration in Eqs. (|22|l 
and H25(l to be the Schwinger-Keldysh contour W along the real-time axis or closed time path (CTP) contour. 
The Schwinger-Keldysh formalism is a powerful method for deriving real and causal evolution equations for the 
expectation values of quantum operators for nonequilibrium fields. The basic idea of the CTP formalism relies 
on the fact that a diagonal matrix element of a system at a given time, t = can be expressed as a product 
of transition matrix elements from t = to t' and the time-reverse (complex conjugate) matrix element from 
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t' to by inserting a complete set of states into this matrix element at the later time t' . Since each term 
in the product is a transition matrix element of the usual or time reversed kind, the standard path integral 
representation for each one can be introduced. However, to get the generating functional we seek, we have 
to require that the forward time evolution takes place in the presence of a source J+ but the reversed time 
evolution takes place in the presence of a different source J~ , otherwise all the dependence on the source drops 
out. 

The doubling of sources, the fields and integration contours suggest introducing a 2 x 2 matrix notation. 
This notation has been discussed extensively in the literature (see [SI d)- However we are going to follow 
Refs. 15 and l_6j and introduce the CTP formalism in our equation of motion by using the composition rule 
for transition amplitudes along the time contour in the complex plane. This way is cleaner notationally and has 
the advantage that all the functional formalism of the previous section may be taken with the only difference 
of path ordering according to the complex time contour Cctp in the time integrations. 

The two-point functions are decomposed as 

G^^M = e,tpMG^^>{t,t') + e,tp{t\t)Gp{t,t'), (26) 

where 



hG1^>{t,t') = {^nt)^W))^ (27) 
hGt^<{t,t') = {^'l{t')v';{t)) , (28) 

with being the fluctuation field and 9ctp{t — t') being the CTP complex contour ordered theta function defined 

by 

e(t, t') for t and t' both on C+ 



, _ I 0{t',t) for t and t' both on G , , 

yctp(t,t)-\ ^^^^ ^ ^, ^+ (29) 

for t on C+ and f on 



With these definitions the matrix indices are not required. When integrating over the second half G , we 
have to multiply by a negative sign to take into account the opposite direction of integration. 

To show explicitly that the prescription for the CTP integration explained above does lead to a well-posed 
initial value problem with causal equations, let us explicitly consider the integral in Eq. H25|) . The integrand 
has the CTP ordered form 

s(t, i")G(t", t') = e,tp{t, t")^?ctp(t", m> (t, t")G> it", t') + e,tp{t, t")e,tp{t', t")s> (t, t")G< {t", t') (30) 
ectp{t\ t)e,tp{t!', t')^< [t, t")G> it", t') + e,tp{t", t)e,tp{t', t")s< {t, t")G< [t", t'), 

where we have omitted the indices because they are not relevant for the discussion. Using the rule for CTP 
contour integration we get 



dt"^{t,t")G{t" ,t') = / dt" {e{t",t')^>{t,t")G>{t",t') + e{t",t")i:>{t,t')G<{t",t')) 

Jo 

+ I dt" {e{t", i')s< {t, t")G>{t", t') + e{t', t")s< [t, t")G<{t", t')) 

dt"^<{t,t")G>{t",t'). (31) 
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li t > t', we have 



/ dt"J:{t,t")G{t",t')^ f dt"{^>(t,t")-Y.<{t,t"))G>{t",t')^ f dt"Y.> {t,t"){G> (t" ,t')~G<(t" ,t')). (32) 
J Jo Jo 

On the other hand, ii t < t' 

dt"^{t,t")G{t" ,t') = f dt"{T?-{t,t")-T.<{t,t"))G<{t\t')- f dt"T.> {t,t"){G> {t" ,t')-G<{t" ,t')). (33) 
Jo Jo 

The above equations are exphcitly causal. 

It is convenient to express the evolution equations in terms of two independent two-point functions which 
can be associated to the expectation values of the commutator and the anti-commutator of the fields. We define, 
following Ref. 16 the functions 

G(f'^''(t,t') = \{Gt^>{t,t') + G'i^<{t,t')) (34) 
G^^'>^\t,t') = ^{G'i^>{t,t')-Gt^<^t,t')), (35) 

where the (F) functions are usually called statistical propagators and the (p), spectral functions. (See Ref. 
|33j ) . With these definitions Ea. (|25|l can be rewritten as: 

U 



KhdtGiP^'it, t') ^ -J [Gi^lfit, t') + GZ\f{t, t')) + (<j>,Mm)GZ^''\t. t')) + (36) 







\tl^f,^^\t^t'')Gif{t'',t'), 



Knd.G^f^'it, t') = J (G^^,^{t, t') + G^^^it, t')) + ^ (c^^^mit))G\f\t, t') + (37) 



Af 

In particular, we define the normal, p, and anomalous, m, spectral and statistical functions as 

Gf''\t,t') ^ pif (t,t') = i((/.l(i)^,(i') + ^,(O^IW), (38) 
Gf'\t,^) ^ p^\t,t')^^{^l{t)^^it')-^^it')^lit))^ (39) 

Glf^\t,t') ^ m(f (t,<')-^(^.(0^,(i') + ^,(i')^,(i)), (40) 

With these relations in place, we now proceed to derive the time evolution equations for the mean field 
and the two-point functions from the CTP 2PI effective action for the Bose-Hubbard Model under the three 
approximations described before. 
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6 HFB Approximation 



As remarked in Section 2 the first order mean field approximation leads to a DNLSE which includes only the 
contribution from the condensate. The HFB equations go beyond the first order approximation and include the 
leading order contribution of They describe the coupled dynamics of condensate and non-condensate atoms 
which arise from the most important scattering processes which are direct, exchange and pair excitations. The 
basic damping mechanism present in the HFB approximation are Landau and Beliaev damping associated with 
the decay of an elementary excitation into a pair of excitations in the presence of condensate atoms ( [251 12tj| l. 
However, these kinds of damping ^ found in the HFB approximation (due to phase mixing, as in the Vlasov 
equation |27j) are different in nature from the collisional dissipation (as in the Boltzmann equation) responsible 
for thermalization processes. Multiple scattering processes are neglected in this approximation. We expect the 
HFB equations to give a good description of the dynamics in the coUisionless regime when interparticle collisions 
play a minor role. 

The leading order contribution of r2 is represented by the double-bubble diagram. Its contribution to T2 
is (f> independent and has an analytic expression of the form 

4'^[G] = dU (G4(t,t)G,?,(i,i) + 2G,,,,(i,i)Q'(i,i)) , (42) 

the factor of two arises because the direct and exchange terms are identical. 

Using the first order expression for r2 in Eqs. H22|l and (|25|l yields the following equations of motion. 



Khdtclylit) = 3^^^, (43) 

EE -j(0^+i(i) + ^ {<l>^Am^{t)+G,'^,{t,t)) c/^Ut) + ^ (^'(OG,,,, (<, i)) , 



KdtG^!; {t,t') = (44) 

^ -J{Gt^,,{t,t') + Gl\,{t.t')) + ^{<l^umi{t)+Gi^{t,t))G'i^{t,t') + 
211 

— {^imc{t)G'^' it, t') + G,?,(t, t)Gf^{t, t')) + zrt,^c(i - t'). 

In terms of the spectral and statistical functions, Eqs. to ||1T|) . and setting N =2, the above equations 
take the form 



ifidtUt) = -J(0,+i(O + '/'.-i(i)) + t^(l0.(t)l' + 2pi''Ht,t))</>,(O + Mr(i,i)'/':(^)> (45) 

-^n^/^\t,t') = U^(t)pf^{t,t') + M.Ut)mf^{t,t'), (46) 

-^n^^p\f{t,t') = Uk{t)p[']{tX) + y^UtW^]{t,t'), (47) 

^h^m^^\t,t') = L,,.(t)mif (t,t') + M,,.(t)p(f (t,t'), (48) 

ih^m^^{t,t') = Uk{t)m%\t,t') + m,,{t)pi'^{t,t'\ (49) 



^we make a distinction between the meaning of the words 'damping' and 'dissipation', the former referring simply to the 
phenomenological decay of some function, the latter with theoretical meaning, e.g., in the Boltzmann sense 
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with 



My(t) 



-J{h+ij + h-i,) + 2(75,, {W)? + f^(t, t)) , 



(50) 
(51) 



The time dependent HFB equations are a closed set of self-consistent equations that describe the coupled 
dynamics of the condensate and non-condensate components of a Bose gas. It can be checked that they preserve 
important conservation laws such as the number of particles and energy. The conservation properties of the 
HFB equations can also be understood by the fact that these equations can also be derived using Gaussian 
variational methods '24'. These methods always yield a classical Hamiltonian dynamics which guarantees 
probability conservation. Because they are local in time they can be decoupled by a mode decomposition. (See 
Appendix 1 for details). 

7 Second order expansion 

7.1 Equations of Motion 
7.1.1 Full second order 

The second order contribution to r2 is described in terms of the setting-sun Fig. (2b) and the basket-ball Fig. 
(2c) diagrams. The basket-ball diagram is independent of the mean- field and is constructed with only quartic 

(2) 

vertices. The setting-sun diagram depends on and contains only three-point vertices. The second order r2 
effective action can be written as 




To simplify the notation, let us introduce the following definitions |16| : 



%(t,t') = --G,,ab{t,t')G'}^{t,t'), 

S„ab(i,t') = -D{t,t')G,.^,{t,t'), 



(53) 
(54) 
(55) 

(56) 
(57) 
(58) 



Dit,t') = cP,,{t)(^^,it')Gtf {t,t')-n,^{t,t'), 



^ijait^t') — ^Gij{t,t')Gij^^{t,t'), 

A,5-„(i,t') = -G,5^(t,i')G,,,,(i,i'), 



With the above definitions we find from Eqs. (|22|l and (|25|) the following equations of motion: 
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+'{w)' J dt''i^^tditX) + ^^^{tX)^Jt)<t>,An)Gt'; {t",t'), (60) 

where S^rFS ^^'^ ^hfb ^^"^ defined in Eqs. 144|l and 145|l . For explicit expressions in terms of p^^'P^ and 
m^F'P^ see Appendix 2. 

7.1.2 2PI-1 /TV expansion 

The 2PI effective action is a singlet under 0{Af) rotations. It can be shown that all graphs contained in an 
0(7V) expansion can be built from the irreducible invariants |16|: (j)^, tr{G") and tr{(l)4>G'^) , with n < J\f. The 
factors of A/" in a single graph contributing to the same 1/Af expansion have then two origins: a factor of J\f 
from each irreducible invariant and a factor of 1 /Af from each vertex. The leading order large Af approximation 
scales proportional to J\f, the next to leading order (NLO) contributions are of order 1 and so on. At leading 
order only the first term of Eq. (|42(l contributes. At the next to leading order level, if we truncate up to second 
order in the coupling strength, the double-bubble is totally included but only certain parts of the setting-sun 
and basket-ball diagrams are: the first term in both of the integrals of Eq. H5,S|I . 

dt,dt,<j,^,{t)<j,^,,{t') [Gtf{t,t')G,,ad'{t,t')Gff{t,t')) + 

f dt,dt'^ (g,,,,, [t, t')G'^f (t, i')G,,d,, (<, t')Gff {t, t')) . (61) 

The equations of motion under this approximation are the ones obtained for the full second order expansion 
but with A = A = 0, and 6 = S. 

In Appendix 2 we explicitly write the equations of motion in terms of the spectral and statistical functions. 
We end this section by emphasizing that the only approximation introduced in the derivation of the equations 
of motion presented here is the truncation up to second order in the interaction strength. These equations 
depict the nonlinear and non-Mar kovian quantum dynamics, which we consider as the primary distinguishing 
features of this work. It supersedes what the second order kinetic theories currently presented can do, their going 
beyond the HFB approximation notwithstanding. For example Ref. presents a kinetic theory approach that 
includes binary interactions to second order in the interaction potential but uses the Markovian approximation. 
In Ref. j35 _ the authors gave a non-Markovian generalization to the quantum kinetic theory derived by Walser 
et. a/.|36| by including memory effects. However in that work symmetry breaking fields, (/) and anomalous 
fluctuations m are neglected. 

7.2 Conservation Laws 

For a closed (isolated) system the mean total number of particles N and energy are conserved quantities as they 
are the constants of motion for the dynamical equations. 

Particle number conservation is a consequence of the invariance of the Hamiltonian under a global phase 
change. The mean total number of particles is given by 
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(n) = T.{^l^^)-T.(\^f+P^'-l) (62) 
= N. 

The kinetic equation for N is then 

J^{m) - pRe[M4^m))+^^^^,{§-,pP (63) 
= 0. 

All three approximations we have considered, namely, HFB, l/J\f expansion and full second order expansion, 
conserve particle number. This can be shown by plugging in the kinetic equation of ^iV(t)^ the equation of 
motion for the mean field and the normal statistical propagator and cancelling terms. It is important to 
note that even though total population is always conserved there is always a transfer of population between 
condensate and non condensate atoms. 

While number conservation can be proved explicitly, to prove total energy conservation is not obvious as the 
Hamiltonian cannot be represented as a linear combination of the relevant operators. It is clear that the exact 
solution of a closed system is unitary in time and hence disallows any dissipation. However, the introduction 
of approximation schemes that truncate the infinite hierarchy of correlation functions at some finite order with 
causal boundary conditions may introduce dissipation j^. 

To discuss energy conservation we can use the phi-derivable criteria |37| which states that nonequilibrium 
approximations in which the self energy S is of the form S^/SG, with $ a functional of G, conserve particle 
number, energy and momentum. All the approximations we consider in this paper are phi derivable and thus they 
obey energy, particle number and momentum conservation laws. For HFB, $ = r2^', for the full second order 
expansion, $ = r2^-' + Fj'^-' and for the second order next to leading order 1/Af expansion, <I> = r2^'' + T^'^^^^ . 
See Eqs. (|24(l . I|42() . (|53|l and H61|l . For a detailed discussion of the complete next to leading order 1/Af expansion 
see Refs. and references therein. 

7.3 Zero mode fluctuactions 

The spectrum of fluctuations above the condensate includes a zero mode. This mode is the Goldstone boson 
associated with the breaking of global phase invariance by the condensate. It is analogous to the collective 
modes which arise in the spectrum of fluctuations around a bubble |38) . The zero mode is essentially non 
perturbative. In linearized theory, it introduces an artificial infrared divergence in low dimensional models. For 
this reason linearized theory is actually improved if the contribution from this mode is neglected all together 
|39| . A different way to deal with the zero mode has been proposed by Gardiner 0D] and Morgan Here the 
theory is written in terms of operators which exchange particles between zero and nonzero modes, conserving 
the total particle number, and one further operator which changes total particle number. The contribution from 
the zero mode is then subtracted by expressing the normal and anomalous densities in terms of the former alone. 
However, from a physical point of view the zero mode exists and is quantum in nature. We may think of it as 
the limit of Gardiner and ZoUer's 0J| "condensate band "when the width of the band shrinks to zero. There 
are both fundamental and practical reasons why isolating and subtracting the zero mode is not as compelling 
in our case as in the problems discussed by Gardiner and Morgan. In the problem we discuss, the initial 
state is a coherent state rather than a proper state of the total particle number. As the total particle number 
is not very high quantum fluctuations in the total particle number are real, and non-negligible. Discarding 
these fluctuations would spoil the integrity of the formalism. Also, because the 2PI formalism goes beyond the 
linearized approximation, the zero mode does not have the impact it has in the linearized formalism and it is 
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not clear that subtracting it necessarily leads to a better approximation. Therefore, in this paper we shall not 
attempt to isolate the contributions from the zero mode. A full non-perturbative treatment in the future is 
certainly desirable. 

8 Numerical Implementation 

8.1 Exact Solution 

The fully quantal solution was found by evolving in time the initial state with the Bose-Hubbard Hamiltonian 
given by Eq. CD,so that \(p{t)) = e-^^*|(/5(0)) with \<f{0)) = e-^'/2e^/^*^ |0)oni5^o To do the numerical 
calculations we partitioned the Hilbert space in subspaces with fixed number of atoms and propagated indepen- 
dently the projections of the initial state on the respective subspaces. A subspace with iV„ number of atoms 
and / number of wells is spanned by states. This procedure could be done because the Hamiltonian 

commutes with the number operator and thus during the dynamics the different subspaces never get 

mixed. The number of subspaces used for the numerical evolution were such that no change in plots of the 
dynamical observables was detected by adding another subspace. Generally for N atoms in the initial state, 
this condition was achieved by including the subspaces between iV — i^/N and N + A^/N atoms. 

8.2 Numerical Algorithm for the approximated solution 

The time evolution equations obtained in Sec. 7 are nonlinear integro-differential equations. Though the 
equations are very complicated, they can be solved on a computer. The important point to note is that all 
equations are casual in time, and all quantities at some later time tf can be obtained by integration over the 
explicitly known functions for times t <tf. 

For the numerical solution we employed a time discretization t — nat, t' — mat, and took the advantage that, 
due to the presence of the lattice the spatial dimension is discrete (indices i and j). The discretized equations 

for the time evolution of the matrices p\^^m ' '^i jnm ^'^'^ 4>in advance time step wise in the n-direction for fixed 
TO. Due to the symmetries of the matrices only half of the (n, to) matrices have to be computed and the values 
Pijnn ~ "^ijnn ~ ^ are fixcd for all time due to the bosonic commutation relations. As initial conditions one 

specifies p]^oo ^ "^Ifoo^ and 0,o. 

To ensure that the discretized equations retain the conservation properties present in the continuous ones 
one has to be very careful in the evolution of the diagonal terms of and take the limit to — > n in a proper 
way: 

\yijn+ln rijnn J \ l^jin+ln 1^ jinn J ' \^^) 
y'Hjn+ln ''Hjnn J ^ y'^jin+ln '''■jinn J ' I'^'^J 

with the positive sign for the statistical propagators, {F)'s, and negative for the spectral ones, {pYs. We use 
the fourth order Runge-Kutta algorithm to propagate the local part of the equations and a regular one step 
Euler method to iterate the non local parts. For the integrals we use the standard trapezoidal rule. Starting 

with n — 1, for the time step n + 1 one computes successively all entries with m = , n, n + 1 from known 

functions evaluated at previous times. 

The time step at was chosen small enough so that convergence was observed, that is, further decreasing 
it did not change the results. The greater the parameter UN/ J the smaller the time step required. The 
main numerical limitation of the 2PI approximation is set by the time integrals, which make the numerical 
calculations time and memory consuming. However, within a typical numerical precision it was typically not 
necessary to keep all the past of the two point functions in the memory. A characteristic time, after which the 



(F,p) _ (F,p) 

(Rp) (F.p) 
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influence of the early time in the late time behavior is given by the inverse damping rate. This time is described 
by the exponential damping of the two-point correlator at time t with the initial time'16'. In our numerics we 
extended the length of the employed time interval until the results did not depend on it. In general, it was less 
than the inverse damping rate. We used for the calculations a single PII 400 MH workstation with 260 Mb of 
memory. For a typical run 1-2 days of computational time were required. 

8.3 Initial conditions and parameters 

To model the patterned loading the initial conditions assumed for the numerical solutions were 4>^q = 
N5io, PijOQ = ^Sij, P^jIq — —iSij and to'^q = fn^j^Q — 0. They correspond to an initial coherent state 
with TV atoms in the initial populated well. 

To study the kinetic energy dominated regime we chose for the simulations three different sets of parameters: 
The first set was chosen to be in the very weak interacting regime, / = 3, iV=:6,J=l and U/J = 1/30. With 
this choice we wanted to show the validity of a mean field approach to describe this regime and the corrections 
introduced by the higher order approximations. The second set of parameters were I = 3,N = 8,J=1 and 
U/J ~ 1/3. In this regime the kinetic energy is big enough to allow tunneling but the effect of the interactions 
are crucial in the dynamics. Comparing with the exact solution we could show the break-down of the mean 
field approximation. 

At the mean field level (using the DNLSE) for a given number of wells, the only relevant parameter for 
describing the dynamics of the system is the ratio UN/ J. For a fixed UN/ J the mean field dynamics is 
independent of the number of atoms in consideration. This is not the case in the exact solution where both 
UN/ J and N are important. As N is increased the bigger the population in the initial coherent matter field 
and therefore we expect a better agreement of the truncated theories with the exact solution. To study the 
dependence of the dynamics on the total number of atoms, the third set of parameters in our solutions were 
chosen to be / = 2, J = 1/2 and fixed NU / J = 4 but we changed the number of atoms from 20 to 80. To 
increase the number of atoms in the calculations we had to reduce the number of wells to two due to the fact 
that the dimension of the Hilbert space scales very badly with N and /. 

9 Results and Discussions 

In Figs. 3 to 7 we show our numerical results. We focus our attention on the evolution of the condensate 
population per well, the total atomic population per well, -I- plf\t,t) — i, the depletion per 

well or atoms out of the condensate, p^P {t, t) — i, and the total condensate population,^- The total 

population is also explicitly shown in the figures to emphasize number conservation. 

The quasi-momentum distribution of the atoms released from the lattice is important because it is one of the 
most easily accessible quantities from an experiment. The quasi-momentum distribution function is defined 
as 

n,(0 = j^e"^-(-^")(<i>IW<i>,W), (66) 

where the quasi-momentum k can assume discrete values which are integral multiples of j^, with / the total 
number of lattice sites and a the lattice spacing. 

The basic features of the plots can be summarized as follows: 

In the very weak interacting regime (Fig. 3) the dynamics of the atomic population per well is almost 
like Rabi oscillations. Notice that even though there are three wells, periodic boundary conditions enforce equal 
evolution of the initial empty ones. In this regime damping effects remains very small for the time depicted 
in the plots. The numerical simulations show a general agreement between the different approaches with the 
exact solution. The effect of including higher order terms in the equations of motion introduce small corrections 
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which improve the agreement with the exact dynamics. This shows up in the plots of the condensate population 
and depletion, where the small differences can be appreciated better. The second order I/JV expansion gives 
an improvement over the HFB and the complete second order perturbative expansion almost matches perfectly 
with the exact solution. In the chiration depicted in the plots of Fig. 3 the total condensate constitutes an 
important fraction of the total population. Regarding the quasi-momentum distribution we observe that similar 
to the spatial distribution where the initial configuration and periodic boundary conditions reduce the three 
well system to a double well one, they enforce equal evolution of the ±^ quasimomentum intensities. The 
A; = and ±^ intensities oscillate with the same frequency as the atomic population per well, both are also 
well described by the approximations in consideration. 

In the intermediate regime we can see the effect of the interactions in the dynamics. They modulate the 

oscillations in the population per well and scatter the atoms out of the condensate. 

1. In Figs. 4 we plot the numerical solution for the parameters I = 3,N = 8,J = 1 and U/J= 1/3. Different 
from the very weak interacting regime only at the very early times all three approximations are close to 
the exact solution. Even though any of thorn arc good after the first oscillation the HFB approximation is 
the only one that fails to capture the exponential decrease of the condensate population. This is expected, 
because even though this approximation goes beyond mean field theory and takes into account the most 
important scattering effects, it includes the effects of collisions only indirectly through energy shifts, and 
breaks down outside the coUisionless regime where multiple-scattering effects are important. In contrast, 
the exponential decay of the condensate is present in the second order approximations. Non local parts 
of the self-energy included in them encode scattering effects responsible for damping. It is important to 
point out that, even though we observe the collapse of the condensate population, the total population 
is always conserved: As the condensate population decreases, the number of atoms out of the condensate 
increases. 

2. Comparing the two second order approaches we observe that the full second order expansion gives a better 
description of the dynamics than the 1/J\f solution only in the regime when the perturbative solutions 
are close to the exact dynamics. As soon as the third order terms start to be important the large 1/Af 
expansion gives a better qualitative description. This behavior is going to be appreciated better in figs. 5 
to 7 as the number of atoms is increased (see discussion bellow). 

We observe as a general issue in this regime that, regardless of the fact that the second order solutions 
capture well the damping effects, as soon as the condensate population decreases to a small percentage 
of the total population, they depart from the exact dynamics: the second order approaches predict faster 
damping rates. The overdamping is more severe in the dynamics of the population per well than in the 
condensate dynamics. The failure can be understood under the following lines of reasoning. At zero 
temperature condensate atoms represent the most "classical" form of a matter wave. When they decay, 
the role of quantum correlations become more important. At this point the higher order terms neglected 
in the second order approximations are the ones that lead the dynamical behavior. Thus, to have a 
more accurate description of the dynamics after the coherent matter field has decayed one needs a better 
treatment of correlations. 

Damping effects are also quite noticeable in the quantum evolution of the quasi-momentum intensities. 
Similar to what happens to the spatial observables, the HFB approximation fails completely to capture 
the damping effects present in the evolution of the Fourier intensities whereas the second order approaches 
overestimate them. 

3. In Figs. 5 to 7 we explore the effect of the total number of atoms on the dynamics. In the plots we show 
the numerical solutions found for a double well system with fix ratio UN/ J = 4 and three different values 

of N: N=20,40 and 80. Wo present the results obtained for the evolution of the atomic population per 
well in Fig. 5, the condensate population per well and total condensate population in Fig. 6 and the 
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quasi-momentum intensities in Fig. 7. To make the comparisons easier we scaled the numerical results 
obtained for the three different values of N by dividing them by the total number of atoms. In this way 
for all the cases we start with an atomic population of magnitude one in the initial populated well. 

In the exact dynamics we see that as the number of atoms is increased the damping effects occur at slower 
rates. This feature can be noticed in the quantum dynamics of all of the observables depicted in the plots 
5 to 7. The decrease of the damping rates as the number of atoms is increased is not surprising because 
by changing the number of atoms we affect the quantum coherence properties of the system. As shown in 
Ref. [3n| , the collapse time of the condensate population is approximately given by tcoU ~ where a is 
the variance of the initial atomic distribution and trev is the revival time which depends on the detailed 
spectrum for the hamiltonian. In the kinetic energy dominated regime trev ~ hN/J (see Ref.|19p. thus 
tcoU For our initial conditions the variance is proportional to cr = so tcoU ~ V^- Besides 

damping rates, the qualitative behavior of the exact quantum dynamics is not affected very much as the 
number of atoms is increased for a fixed UN/ J. 

The improvement of the 2PI approximations as N is increased, as a result of the increase in the initial 
number of coherent atoms is in fact observed in the plots. Even though the problem of underdamping in 
the HFB approximation and the overdamping in the second order approaches are not cured, as the number 
of atoms is increased, we do observe a better matching with the full quantal solution. The 1/Af expansion 
shows the fastest convergence. Perhaps this issue can be more easily observed in the quasi-momentum 
distribution plots. Fig. 7. The better agreement of the 1/J\f expansion relies on the fact that even though 
the number of fields is only two in our calculations the 1/M expansion is an expansion about a strong 
quasiclassical field configuration. 

10 Conclusions 

In this work we have used the CTP functional formalism for 2PI Green's functions to describe the nonequi- 
librium dynamics of a condensate loaded in an optical lattice on every third lattice sites. We have carried out 
the analysis up to second order in the interaction strength. This approximation is introduced so as to make 
the numerical solution manageable, but it is sufficient to account for dissipative effects due to multiparticle 
scattering that are crucial even at early times. Our formulation is capable of capturing the salient features 
of the system dynamics in the regime under consideration, such as the decay of the condensate population 
and the damping of the oscillations of the quasi-momentum and population per well unaccounted for in the 
HFB approximation. However, at the point where an important fraction of the condensate population has been 
scattered out, the second order approximations used here predict an overdamped dynamics. To improve on 
this a better treatment of higher correlations is required. One might try to include the full next to leading 
order large Af expansion without the truncation to second order as done in Ref. (16) but it is not obvious that 
this will lead to the required improvement. Alternatively, one may try to adopt a stochastic approach, but the 
challenge will be shifted to the derivation of a noise term (which is likely to be both colored and multiplicative) 
which contains the effects of these higher correlations and the solution of the stochastic equations. We hope to 
address this aspect of the problem in a a future work. 

Even though, as is clear in this paper, the second order 2PI approximations fail to capture the fully correlated 
dynamics in the system, it has been proved to work at intermediate times when correlations are not negligible 
and standard mean field techniques fail poorly. Because of its success in describing moderately correlated 
regimes, the second order approximations could become a useful tool for describing experimental situations as 
the collapsing or colliding condensates experiments where striking dynamical behavior such as collisional loss 
of condensate atoms has been observed. 

In summary, we have presented a new approach for the description of the nonequilibrium dynamics of a 
Bose-Einstein condensate and fiuctuations in a closed quantum field system. The formalism allows one to go 
beyond the well known HFB approximation and to incorporate the nonlinear and non-Markovian aspects of 
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the quantum dynamics as manifest in the dissipation and fluctuations phenomena. The 2PI effective action 
formahsm provides a useful framework where the mean field and the correlation functions are treated on the same 
footing self-consistently and that respects conservations of particle number and energy. The CTP formalism 
ensures that the dynamical equations of motion are also causal. In their current form scattering terms non 
local in time are hard to estimate analytically and numerically demanding. However, this systematic approach 
can be used as a quantitative means to obtain solutions in different regimes and make comparisons with kinetic 
theory results where a Markovian approximation is assumed. 
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Appendix 1: Mode Expansion of the HFB Equations 

To decouple the HFB equations we apply the well known Bogoliubov transformation to the fluctuation field: 

Vj{t)=J2ut{t)a,-v*\t)al, (A-1) 

where (aq,a]j) are time independent creation and annihilation quasiparticle operators and all the time de- 
pendence is absorbed in the ampHtudes {u'^ (t) , v*'' (t)} . To ensure that the quasi-particle transformation is 
canonical, the amplitudes {uj{t),v*'^{t)} have to fulfill the following relations^: 



^t^KO^fW-^KiX'W = Sik, (A-2) 

i 

J2ul{t)v^{t)-vUt)u^{t) = 0. (A-3) 



In the zero temperature limit, where ia^^ak) = 0, the statistical and spectral functions take the form 





q 


+u]{t')unt)), 




q 






q 


+ u]{t')v*%t)), 









q 



(A-4) 
(A-5) 

(A-6) 
(A-7) 



^The HFB approximation can be seen as a quadratic approximation of the Hamiltonian in which third and quartic terms are 
reduced to linear and quadratic forms by factorizing them in a self-consistent Gaussian approximation. Any quadratic Hamilto- 
nian can be diagonalized exactly and this is a<;hieved by transforming to a quasi-particle basis, {&q}. The requirement of the 
transformation to be canonical means that it preserves commutation relation and leads to bosonic quasi-particles. 
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Replacing Eqs. HA-4|I - HA-7|I into Eqs. and using the constraints (|A-2|) - (|A-3|I we recover the standard 

time dependent equations for the quasiparticle amphtudes 24 : 

ihdtm = -J(0,+i(t) + + U{\m? + 2p^p{t,t))cj>S) + Un4P{t,t)^*{t), (A-8) 

= -Jiu'!^,{t)+uU{t)) + 2Umt)\'+P^\t,t))u!{t) - U{m*f \t,t) + ^:{tf )v'l{t){A-%) 

-ifi^^vKt) = ~j{vt^,{t)+vum + 2u{\m? + P^Pit^tMit) - u{mPHt,t)+m^)<{tik-iQ) 

Equations (|A-8(I - l|A-10(l correspond to a set of 7(2/ + 1) coupled ordinary differential equations, where 
/ is the total number of lattice sites. They can be solved using standard time propagation algorithms. Once 
the time dependent quasiparticle amplitudes is calculated we can derive the dynamics of physical observables 
constructed from them as a function of time, such as the average number of particles in a well ni(t), etc. 

Appendix 2: Second order equations of Motion 

Hereunder we explicitly write the equations of motion of the 1/M and full second order approximations. 
To simplify the notation let us introduce the functions 

n[P[i,g] = f(f)(i„t,)glf (t„i,)-i(f(f)(t„t,)g(f(i„t,)), (B-i) 

n[f[i,g] = f^)(t„t,)gg'(t„t,)+f<f)(t,t')g!f (B-2) 

Using the spectral and statistical functions and setting J\f equal to 2 the equations of motion derived in Sec. 7 
can be written as: 

1. Full second order expansion 

If we take the complete contribution of the setting-sun and basket-ball diagrams the equations of motion 
take the form 



ihdtA, = -J(0,+i(t,) + '^.-i(t.)) + C/(l0.P + 2pifV, + C^4fV* (B-3) 



-2C/2 J2 f^dtk [q^.niZ^ [p, p*] + 4>k^^:^ [m, m*] + <t>in^P^ [m, p]) p^' 

A; ''^ 

-2C/2 ^ USI^ [p, p*] + 0*17^^) [m, m*] + 0,17^^) [m\ p*]) m^f 

+2C/2 dtk [cj^^nfP [p, p*] + ) [m, m*] + (j^l^^^ [m, p]) 

k 

+2C/2 E r [^l^P [/^' + ^l^T + 4>k^T Vn\p*]) ^ 



(p) 

ki ' 
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kj 



-Jip^Uu,tj) + p^Uu, tj)) + 2C/(|</,,|2 + p(f Vlf ^ + U{m*P + <l>f)m^^ (B-4) 

-2[/2 ^ rfi, (0(P) [p, A] + 2</.*<^, [p, p*] + l^l,^) [m, m*])) p^f 

-2C/2 5^ dtu {2<l^icl^l^\'^ [m*,p] + <j>,<l>,n^^ [m*,m*] + 2<j>*<p,n^^ [p*,m*]) m 

k 

-2f^' E {^"fk K, T] + 0*0^ (20(,^) [p, p*] + 20(,^' [m, m*])) mg 

+2C72 ^ dtu {<l^i4>l^fu^ [P, P] + 2<^,<^fel^ir [p, m*] + 2<^*<^*^l(f [m, p]) p^J 

J2 r dtu (Olf ) [p, A] + 2</.*0, (17(f) [p, p*] + l^^P [m, m*]) ) pg^ 
fe 

fc -^0 ^ 

+2C/2 5] rftfe (n(f ) [m*, T] + 20*^;^ ) [p, p*] + 0(f) [m, m*])) m^f 



"fei 



+ P\tj{ti^ tj)) + 2Um^ + p(f Vlf + U{m:P + ,^)n^ (B-5) 



-2C/2 ^ / dtfe U,^l^^^ [p, p] + 2</,,</,,0(^) [p, m*] + 2</,*,^^0(^) [m, p] 

-2C/2 t dtu (fiL'^ [P, A] + 2<p*cPu (f^L"' [P> Pi + f^L"^ [m, m*])) pg 
fe -^'j 

-2C/2 ^ /■*' (2</),</)^0(^) [m*, p] + </),0,o(^) [m*, m*] + 2</,*,^,0(^) [p*, m*]) mg 
fe 

-2f/2^ Tdt, (n^^^ [m*,T]+2m (nl^^ [p,P*]+n^^^ Km*])) mg, 
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-J{m^l{U,W + mZUtutj)) + 2Um^ + pif )mlf + [/(mlf + </)?)plf (B-6) 
-2U' ^ /*" dU [m, p*] + m^ll^ [m, m] + 2^,^ln[^'> [p, m]) p^f 



-2U' / rftfc U*<l^k^^^ [p^p*] + 2^^^^ [p^m] + 2<p,cl^,n^^ [m*,p*]) m 

^ t dtk (fil,^) [p*, A] + 20,0^ [p, p*] + 1^1,^) [TO, TO*])) mif 

+2f/2 ^ r difc (2<^*</),^llf [m, p*\ + <^*<^:i^if [m, to] + 2</,,<^*fi(f [p, m]) pg^ 

+2C/2 ^ dt^ (n<i^ [m, T] + 2</.,0, ) [p, p*] + n<i^ [m, m*]) ) p^J 
fc 

+2t/2 ^ dtk {<P*h^^r [p*,P*] + 2<t>l<t^ldi^ [p*M + -^MMP [m*,P*]) 

k 

+2U^ ^ r dt, ) [p*, A] + 2</,,0^ (O^^) [p, p*] + fi^f [m, m*])) TOg\ 



- JK%-(ti,t,-) + mY\.{U,ti)) + 2[/(|0,|^ + p,\^OKf + U{m\t' + <}>t)p\f (B-7) 
+2C/2 ^ dt, (2</.*</,,fi(^' [m, p*] + ,^*</.^0(^^ [m, m] + 2</,,<^^fi(^' [p, mj) pg 



+2C/25] / dtfe(o(,^)[m,T] + 0,<A,(2f2(,^)[p,p*] + 2fi(,^)[TO,TOl))pg 

fe 

+2U' ^ /*' dtfe (<^*'/'fci^L'^ Pi + 2</'*<^fcOL^) [p*, TO] + 2</,,0,o(,^) [TO*, p*]) TOg 
fc 

+2[/2 ^ r di, (nl^^ [p* , A] + 2<p,rk (f^L'^ [P, Pi + ^L'^ [m, m*]) ) mg) , 



Af'') = S}[f''^[p,p*]+2n[f'''^[m,m*], (B-8) 
T(f^) = 2Q\f'''^[p,p*]+Q\f'''^[m,m*]. (B-9) 
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2. Second order l/M expansion: 



indt,(t>i = -J(</'i+l(^^) + 0._l(^^)) + c/(|</'.r + 2pi^),^, + ^7m^fV 



(B-10) 



k 

-U^V] r dtk 
k -^0 

k Jo 

k -^0 

"Y^ r dtk 
k Jo 

Y dtk 
k Jo 

V / ' dtk 

k Jo 

L. Jo 



+U 



+U 



fe "'o k •'0 

J{P%\j{ti,h) + Pt\j{tu tj)) + 2Ui\ct>f + f^P)pP + U{m*P + </)f )m,f (B-11) 

<t>^<l^l^\^ [P.P] + <t>^<l^k^\^ + 0:0:01,^' [m,p\) 

[p, n] + 4>*4>k {^^^^ [p, p1 + ^'t^ [m, m*]) ) p(f 

K,n] + </,*</,^ [^^J;} [p,p*] + 2n(,^) [m,m*])) mif 

'^i'/'fef^lf ^ [P, P] + '/'^'^fe^ir b, ml + ^r^^f^lf ^ [m, P\) pif 
l^if [p, n] + 0:0, (217(f [p, p*] + ) [m, m*])) pg^ 
</.,</)^n|f K,p] + 0,</.feOif [m\m*] + ^*cf>kn^^ [p*,m*]) 
0(f K,n] + c/-:,^;^ [p,p*] + 20^^ [m,m*])) mip, 



-Apl'^ijitutj) + pfijituW) + 2Um' + f^')f^' + Uimir' + )m^f (B-12) 
-U' E d*'^ {'t>icl>l^fk [P^ P] + 't>i^k^fk [P^ "^1 + 't>Wk^'-£ [m, p]) pg 



kj 



kj 



-U'Y / dtk{n<^Hp,m+<l>:<l>k{'^^'£^[P,P*]+^'£Hrn,m*]))pi 

k J^i 

-u' E f '^^^ (<^i<^^"l^^ K' + <^i'^/^"l^^ K' ™1 + 't>*4>k^^^ [p\m*\) 

k J*3 

-t/^E r^*" {^ik {m*,U]+<p:<l>l (nlf [P,P1 +2^7^^) [m,m*])) 



ip) 



"^kj 



(p) 

kj ' 



k ""-J 
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indum<^^ = -J(m\P,.{tutj)+mZUtutj)) + 2U{\<f>f + p\[^)m^^ (B-13) 
k "'o 

E dt, {n'^^^ [m, n] + 0,<A, [p, p*] + 2n^^ [m, m*])) p^f 

E 'itk [p*, n] + 0,</,^ (21l(,^) [p, p*] + l^l,^) [m, m*])) mif 

+u' E ^^'^ (<^:<^fe^ir Pi + •/'I'/'fc^ir H + 0i<^:^ir ip, h) pI? 

^ •'0 



^ dt, (l^lP [m, n] + 0,</., (217(f [p, p*] + [m, m*]) ) pg^ 
^ dtfc (</.:</.fcl^if [p1 Pi + 0r</-^o(f ) [p*, m] + </-,0,fi(f [m*, p*]) mi^) 



+[7^5^ / 'dtfc [p*,n] + <^,0;^ (20(f [p,p1 + fiir [m,m*]))mg: 



+U' r dtu U*<t>k^'-^ [m, Pi + cl>*cl>lil'i [m, m] + ^../-^fJ^^) [p, m]) p^J 



+U'Y[' dtk {nl^^ [m, n] + (j^l,^) [p, p*] + 2^^L'^ [m, ml)) pi? 
+C/2^ r dt, (<^:<^,l^L'^ [pIpI +<^:<^:^^^^ [p*,m] +<^,0,i7(,^' [m*,p*]) 

k ^^i 

+U^Y. r '^^^ [pin] +<^i</': (2f^L^^ [P,Pl [m,m*])) 



^kj ' 



with 



'^^ = f^(f [p, p*] + J^if ) [m, m*] . (B-15) 
In the above equations we have simphfied the notation replacing (pi-itk) by (pf. and mkj{tk,tj) by nikj- 
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Figure 1: Comparisons between the exact and the DNLSE solutions for six atoms and three wells. The time 
is given in units of fi/ J. Top panel, strongly correlated regime (7 = 12); Middle panel, intermediate regime 
(7 = 2) and bottom panel weak interacting regime(7 = 0.2). The solid line is the DNLSE prediction for the 
population per well: |?!'o(OP ^^^^ 2(^)P i^'^^ 1^1- EIj ^^le triangles are used to represent the exact solution 
for the population per well calculate using the Bose Hubbard Hamiltonian (Eq. (<i>Q$o), ($| 2*^1,2)- The 
pentagons show the condensate population per well calculated from the exact solution: |(<&o)P and |(<&i,2)P- 
Due to the symmetry of the initial periodic conditions the curves for the i = 1 and 2 wells are the same in all 
depicted curves . 
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Figure 2: Two-loop (upper row) and three-loop diagrams (lower row) contributing to the effective action. 
Explicitly, the diagram a is what we call the double-bubble , b the setting-sun and c the basket-ball. 



29 




Figure 3: Comparisons between the exact solution(solid line), the HFB approximation (boxes), the second 
order large N approximation (pentagons) and the full 2PI second order approximation(crosses) for the very 
weak interacting regime. The parameters used were 7 = 3, -/V = 6, J = 1 and U / J ~ 1/30. The time is given 
in units of h/ J . In the plots where the population, condensate and depletion per well arc depicted the top 
curves correspond to the initially populated well solutions and the lower to the initially empty wells. Notice 
the different scale used in the depletion plot. In the momentum distribution plot the upper curve correspond 
to the fc = ±27r/3 intensities and the lower one to the A; = quasi- momentum intensity. 
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Figure 4: Comparisons for the case I — 3, N — 8, J — 1 and U/J= 1/3. The time is given in units of h/J.In 
the plots the abbreviation 1st is used for the initially occupied well and 2nd for the initially empty wells. In the 
quasi- momentum plots 5 = 27r/a is the reciprocal lattice vector with a the lattice spacing. 
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Figure 5: Comparison between the evolution of the atomic population per well for / = 2, J = 1/2, NU / J = 4 
and N = 20, 40 and 80. Time is in units of U/ J . In the plots PI states for the atomic population in the initially 
populated wells and P2 for the population in the initially empty wells. The number of atoms N is explicitly 
shown in each panel. To make easier the comparisons all the plots are scaled by N such that the number of 
atoms is always normalized to one. 



32 




Figure 6: Time evolution of the condensate population per well and the total condensate population. Same 
parcmctcrs than Fig. 5. Time is in units of Ti/J. In the plots CI states for the condensate population in 
the initially populated well, C2 for the condensate in the initially empty one and CT for the total condensate 
population. All the plots are scaled to have the total number of atoms set to one for all the different N cases. 
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Figure 7: Dynamical evolution of the quasi- momentum intensities. The parameters used were 7 = 2, J = 1/2, 

NU / J = 4 and N = 20, 40 and 80. Time is in units of h/ J. In the plots ko denotes the fc = quasi-momcntum 
component and kl the k = Tr/a one (a the lattice spacing). The plots are scaled to set the total quasi-momentum 
intensity to one for all the different N cases. 
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